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New boundary conditions for granular fluids 
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We present experimental evidence, which contradicts the the standard boundary conditions used 
in continuum theories of non-cohesive granular flows for the velocity normal to a boundary u-n = 0, 
where h points into the fluid. We propose and experimentally verify a new boundary condition for 
u-n, based on the observation that the boundary cannot exert a tension force F& on the fluid. 
The new boundary condition is u • ft = if Fb • n > else ft • P • n = 0, where P is the pressure 
tensor. This is the analog of cavitation in ordinary fluids, but due the lack of attractive forces and 
dissipation it occurs frequently in granular flows. 

PACS numbers: 45.70.-n, 51.30.+j, 51.10.+V, 64.70.Hz 



I. INTRODUCTION 

Granular materials, collections of particles which dis- 
sipate energy through inter-particle interactions have 
tremendous technological importance and numerous ap- 
plications to natural systems. They also represent a seri- 
ous challenge to statistical and continuum modeling, due 
to the small number of particles and dissipation. A fun- 
damental understanding of granular systems, comparable 
to our current understanding of fluids and solids, does not 
exist today [T| but would have far reaching implications 
across many industries. 

The idealized granular system, a collection of identical 
inelastic hard spheres, is a laboratory scale analog of the 
canonical system in kinetic theory — the elastic hard 
sphere gas. Using this analogy, an inelastic version of the 
Boltzmann-Enskog equation have be derived 0] as well 
as an inelastic version of continuum equations for mass, 
momentum, and energy balance Q, [0, 0, H, EH, 0, nearly 
identical to the Navier- Stokes equations. These type of 
equations can produce a quantitative description of dilute 
granular flows with an accuracy of 1% [9[. Boundary 
conditions for velocity and granular temperature have 




FIG. 1: photographs of a 2D granular layer under free fall 
conditions (top row) and vibrated (bottom row). The layer 
is initially in contact with the boundary (first column), but 
then moves away (second and third column). 



been derived to complete these equation [TO, EH, EH, [l3[ , 
which assume that the fluid never leaves the boundary. 

Normal fluids remain attached to boundaries due to 
two effects — adhesion and pressure (external and in- 
ternal). External pressure is unimportant and adhesion 
is not present in non-cohesive granular materials leaving 
the internal pressure alone to hold a granular fluid to a 
boundary. However, if the force on the granular fluid ac- 
celerates the fluid to a velocity greater than the average 
root-mean-squared particle velocity (square root of the 
granular temperature), then the boundary and the gran- 
ular fluid will separate. This is the analog of cavitation 
in normal fluids. While cavitation is unusual in normal 
fluids, low granular temperature at boundaries due to in- 
elastic loss combined with the lack of external pressure 
and adhesion makes cavitation common in granular fluids 
(see Fig. [1]). Common examples include rotating drums 
and vibrated layers. Many features of pattern formation 
in vibrated layers [14] can be understood in terms of the 
time that the layer spends off of the plate [15[. It dif- 
ficult to explain this dependence with models which do 
not allow the granular fluid to leave the boundary. 



II. NEW BOUNDARY CONDITION 

The standard boundary condition [TO, EH, E3, EH for 
the velocity normal to the wall is 



(u — v) • h = 0, 



(i) 
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for the fluid velocity u on a boundary with inward unit 
normal n and wall velocity v. In order to enforce this con- 
dition the boundary must exert a force F5 on the gran- 
ular fluid. Because there is no attractive force between 
the grains and the boundary F5 • n > 0. When the force 
needed to maintain zero velocity is negative the bound- 
ary condition must be changed. Since the boundary force 
then must be zero a no stress condition ft • P • n = is the 
logical choice. As the normal wall velocity will no longer 
be zero we must have a another condition to enforce the 
no particle flux condition at the boundary. Analytically 



FIG. 2: Plot showing the lower boundary of the vibrated 
system as a function of time. 94 cycles at 10Hz are overlaid. 
The horizontal line shows the diameter of a particle. 



more work needs to be done to understand this new con- 
dition on the flux = J pu ■ dS. Numerically, we treat 
the density and velocity as centered in the grid and the 
flux is defined on the boundaries of the grid. In that case 
the density and velocity in the cell nearest the boundary 
are both non-zero, but the flux on the wall is zero (see 
(§HVj) for details of the flux condition). Altogether we 
have 

f (u-v)-n = if F 6 -n>0 , , 

\ n • P • n = 0, $ = else F 6 • n < 0. 1 ] 



III. EXPERIMENT 

We have performed two types of experiments to eluci- 
date the role of boundary conditions in continuum equa- 
tions of motion for granular materials — free-fall and 
vibrated. We place TV (51-117) spherical stainless steel 
ball bearings of diameter D = 3.175 mm in a container 
17.5 D wide by 20 D tall by 1 D deep as shown in Fig. 
[TJ We define the number of rows R = N/17, where 17 
is the number of particles to fill an entire row in a close- 
packed structure. We control a thin plunger, which slides 
through a slot in the bottom of the cell. For the free-fall 
experiment [Fig. []Jtop row)] the plunger pushes all of the 
particle to the top of the cell to start the experiment. The 
particles are then vibrated in order to compact them into 
a perfect hexagonal packing [Fig. [TJtop-left)]. Then the 
plunger is pulled downward with an initial acceleration of 
4g, where g is the acceleration of gravity. The particles 
are then free to fall under gravity. This process is un- 
der computer control and can be repeated many times. 
In the vibrated experiments [Fig. [Tfbottom row)] the 
plunger oscillates sinusoidally, but is offset from the bot- 
tom of the container to produce a half-sine-wave excita- 
tion. The measured position of the plate is shown in Fig. 
El 94 oscillations are superimposed to show the repeata- 



FIG. 3: Color online: Steps in the spatial and temporal aver- 
aging process for density field, (top) Modified Voronoi cells, 
(middle) footprint preserving spatial average over 2 particle 
diameters, (bottom) average over phase (time). 



bility of the drive. The excitation is characterized by the 
non-dimensional maximum acceleration T = A(2irf) 2 /g, 
where A is the maximum amplitude of the plunger and 
/ is the frequency. Using high-speed digital photography 
we measure the positions of the plunger and all of the 
particles in the cell with a relative accuracy of 0.2% of 
D or approximately 6/im at a rate of 840 Hz. We track 
the particles from frame to frame and assign a velocity 
to each one, typically ~ D/b per frame. 

To compare with the continuum theories through our 
simulations (see Section llVj) , we must average the parti- 
cle trajectories to created density, velocity, and temper- 
ature fields. We are interested in studying the flow as 
the grains separate from the boundary so we focus on 
the density fields. To create average fields we run each 
flow condition 90-100 times. In each frame, we construct 
a modified Voronoi cell around each particle and fill it 
with a uniform number density of 1 particle divided by 
the area of the cell. We modify the Voronoi cells on the 
border by including only the area that overlaps with a 
sector created by two lines both tangent to the particles 
edge and tangent to two other neighbors less than two 
particle diameters away. This is the union of the con- 
vex hull formed by all all particle pairs whose centers lie 
within 2 particle diameters of each other. The results of 
such a construction are shown in Fig. 09(top). Next, we 
apply a mask preserving spatial average over 2 particle 
diameters [Fig. [^middle)]. The mask is formed by the 
union of all of the modified voronoi cells. Finally, we bin 
the frames according the phase of the cycle and average 
each bin [Fig. [3^bottom)]. From these images we have 
determined that the flow is uniform in the horizontal x 
direction for both free-fall and vibrations. Therefore, x 
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to obtain one dimension density field along the vertical 
z direction. Space-time intensity plots of these densities 
are shown in Fig. [U 

As a further quantitative comparison, we calculated 
the discharge rate for the free-fall case [Fig. Ufa)] and 
the acceleration of the center of mass acM minus g for 
all cases [Fig. E{b-d)]. M(a C M - g) = F 6 , where M 
is the total mass of the fluid. In Fig. [5jb-d) we plot 
a\) = 1/MFfc • z = (acM — g) • z for both the free fall and 
the vibrating case. 



IV. SIMULATION 

We numerically integrate the continuum equations for 
the balance of mass, momentum, and energy have been 
derived for two-dimensional (2D) granular flows [3L Il6j|: 



p( w+ u.Vu) 

dT 

/%+u-W) 







V P-g 

V q-P :E-7, 



(3) 
(4) 
(5) 



where p is the mass density, g is the gravity vector, T 
is the granular temperature, q is the heat flux vector, 
Eij = 7}(diUj+djUi) are the elements of the symmetrized 
velocity gradient tensor E, and 7 is the temperature loss 
rate. The constitutive relations for the pressure tensor P 
and heat flux q are 



and 



P = (P - 2ATrE)I - 2/i(E - (TrE)I) (6) 



q = -«VT, (7) 



where Tr denotes trace and I is the unit tensor. The 2D 
equations close [16| with the equation of state, which is 
the ideal gas equation of state with a term that includes 
dense gas and inelastic effects, 



P = pT[l + (l + e)G(v)}, 



(8) 



is the 



where e is the coefficient of restitution, v = 
volume fraction, a is the diameter of the particles, m 
is the mass of the particles, and G{v) = vg(y,o), where 
g(y, <j) is the value of the radial distribution function at a 
distance of one particle diameter. We use a temperature 
dependent 



e(T) 



l-(l-eo)y/T/n 
eo 



1/5 



ifT < To 
otherwise. 



(9) 



This mimics the experimental evidence that eo goes to 
zero for small velocities. We use a form for G{y) de- 
veloped by Torquato [17[, which is an analytical fit to 
molecular dynamics simulation at high v and the Carna- 
han and Starling [18] geometric series approximation to 



the first few viral coefficients at low v. We have changed 
the functional form of G(y) and the details of the flow 
change, but not the qualitative behavior. The bulk vis- 
cosity, 



A = 2pa\ -G(y), 

7T 



(10) 



the shear viscosity, 



M = 8 L G(y) 
the thermal conductivity, 



2 + (1 + -)£(*,)], (11) 

7T 



3+(^ + ;)GH], 

4 7T 



2 v G{v, 

and the temperature loss rate per unit volume, 
4,C(,)T3/ 2 
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(12) 



(13) 



We use a second-order-space adaptive first-order-time 
finite-differencing scheme to integrate these equation. We 
use off-center-differencing at the boundaries to maintain 
second-order accuracy. The code uses an adaptive time 
step based on a modified Courant condition combined 
with an empirically determined maximum and minimum 
time step. We used a 5x161 point x-z grid, with periodic 
boundary conditions in the x direction. We use the tem- 
perature and tangential velocity boundary condition for 
smooth surfaces found in reference |12j] . We can alternate 
between our new ([2j) and the standard v • ft = bound- 
ary conditions for the normal velocity. To maintain mass 
conservation we use a special differencing procedure for 
the mass continuity equation, which is second-order ac- 
curate for both the new and old boundary conditions and 
is exactly conserving. In one dimension we label the den- 
sity and velocity pi and Wi for i = to N. We then define 
N + 2 fluxes §i = (pi-i + pi)(wi-i + Wi)/4 for i = 1 to 
N, and $0 = $at+i = 0. Then the change in density at 
position Z{ = iAz is Api = (3>i+i — $i)At, where At is 
the current time step. From these definitions it follows 

that Y^iLo = an d J2iLo Pi = constant, regardless 
of the value of w and p at the boundaries. The only 
boundaries in the system are the actual walls. We let 
the free surfaces develop on there own simply as very low 
density regions. As the density of a gas becomes so dilute 
that the mean free path is comparable to the size of the 
container the transport equations must be modified. We 
define a transport cutoff at a density of 



"0 



1 



6I0V2 



5.89x10" 



(14) 



where Iq = 20 is the dimensionless cutoff mean free path. 
At lower densities than this the viscosity and thermal 
conductivity are multiplied by the factor 
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(15) 
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FIG. 4: Color online: Space-time plots of the volume fraction for a) free-fall and b) vibrated. Experiments are shown in the 
first column, simulations with the new boundary conditions in the the second column, and simulations with the old boundary 
conditions in the third column. In a) the upper white line is the top of the cell. The lower line show where the number flux 
is calculated in Fig. [5ja). In b) the line is the position of the bottom plunger. The separation of the fluid from the plate can 
clearly be seen in the experiment and the simulation using the new boundary conditions. 



For v » v>o x is close to unity, but as the density de- 
creases the transport coefficients go smoothly to zero. 
This ad-hoc factor is used only to prevent the code from 
diverging at extremely low densities. In these simula- 
tions, we have varied from 10-100 and the results are 
almost the same except for slight deviations in the very 
low-density regions. This low density regime does not ef- 
fect the main flow as the momentum and energy are also 
nearly zero unless the temperature or velocity diverge, 
which is prevented by this approach. 

We use the following parameter for all of the simula- 
tions presented here: eo = 0.7 and TO = 1. Both of these 
parameter effect the energy loss. If the energy loss is not 
great enough the separation does not occur. These values 
were not formally optimized, but several different value 
of eo between 0.6 and 0.9 were tried, but anything be- 
low 0.8 gave approximately the same flight-time for the 
vibrated case. The temperature and tangential velocity 
boundary conditions require two additional parameters, 
e w = .5 and fi w = 0, the wall coefficients of restitution 
and friction. These parameters only have a small effect 
on the flow. 

For the free-fall case we solve the equation in the lab 
frame. To create an initial condition we set gravity up- 
wards to push the fluid to the top of the container. Then 
we slowly decrease gravity to zero and then abruptly 



change the value of g to —g. The results of the hori- 
zontally averaged density are shown as space-time plots 
in Fig. [Hfor both the new and old boundary conditions. 
For the vibrated case we solve the problem in the refer- 
ence frame of the bottom plate. We determine the accel- 
eration of the plate from the experimentally determined 
position of the plate Fig. [2j The second derivative is 
quite noisy so the signal is convolved with a 5 ms top- 
hat function. The signal is then corrected so that the 
average acceleration is zero. The initial condition is a 
uniform distribution and wait for the density to reach a 
stable periodic state. This typically takes 5-10 oscilla- 
tions. The final cycle is shown in Fig. [H for both types 
of boundary conditions. 



V. RESULTS 

Comparing the first two columns in Fig. [4] (experi- 
ment and new boundary conditions) to the third there 
is a striking effect. In the free-fall experiment the zero 
velocity boundary condition holds the fluid to the top 
surface long after all of the grains in the experiment and 
the new boundary condition simulation have left the do- 
main of interest. To see this in a quantitative way we look 
at the rate the particles pass the lower line in the figure. 
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FIG. 5: Comparison of experiment (connected points) and simulation with new boundary conditions (solid) and old boundary 
conditions (dashed) for free-fall (a,c) and vibration (b,d). Plots of a) Number flux of particles passing the lower horizontal line 
in Fig. [4fa) as a function of time and b-d) acceleration of the center of mass of the granular fluid. 



The result is shown in Fig. [5](a). Initially the two simu- 
lation follow more closely due to the lack of dispersion in 
the experiment. The dispersion in the simulation is not 
due to temperature, but an artifact of the truncation er- 
ror in the mass conservation equations. Shortly ( 20 ms) 
all three curves meet, but by ( 30 ms) the old boundary 
conditions simulation begins to diverge, and from there 
on the behavior is radically different. As the last of the 
material flows from the experiment and new BC simula- 
tion there is another small lag, but the old BC simulation 
continues to flow for another 30 ms. Figure [5^c) shows 
that there is a negative acceleration induced by the old 
BCs. Both the experiment and the new BCs show zero 
acceleration throughout the time of the experiment. 



In Fig. |H[b) the layer is stuck to the boundary for the 
old BCs. There is no clear flight time. From the layer 
acceleration in Fig. [5jb) we see that there is still some 
kind of a collision for the old BCs but it occurs earlier. 
Further, there is a clear region of negative acceleration, as 
seen in the expanded axis of Fig. [Sfd) of nearly half the 
value of gravity. The small negative acceleration in the 
experiment is likely due to friction with the walls. The 
take off time and collision time shown by the two vertical 
line are the same for the experiment and the new BCs. 



VI. CONCLUSIONS 

The new boundary conditions provide a significant im- 
provement in simulating these flows. Separation (cav- 
itation) is a very common situation in granular flows, 
and this type of boundary condition is needed to cap- 
ture the basic features of the flows presented. Without 
this type of boundary condition many qualitative features 
may be missed, and there is no hope quantitative agree- 
ment. There are still a number of open question includ- 
ing: what should happen in corners, and how can these 
BCs be treated analytically? Previously we have shown 
that molecular dynamics (MD) simulation pjl (2(j are 
capable of quantitatively reproducing the flow in dense 
vibrated pattern forming system, and MD and continuum 
simulation can reproduce dilute flows in which granular 
shocks form @. However, this work represents a first, big 
step in getting quantitative agreement using continuum 
simulations in dense vibrated pattern forming systems. 
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